function[w_var_full,w_lp_full,inv_rmse_var,inv_rmse_lp,sum_inv_rmse]=Weight_CVforecast(data,p,testsize,Hstepforecast)
% Weighing based on Cross-validation forecast
% Inputs: data, p (lag-length), testsize, and H-step forecast


% Step 1 Forecast preparation
[T,n]=size(data);
T_half=testsize;       % Max size of test set
H_step_forecast_max=Hstepforecast; % Forecast horizons

no_sample_of_fullforecasts=T_half-H_step_forecast_max+1; % Number of samples having full forecasts
Total_forecast=H_step_forecast_max*(H_step_forecast_max-1)/2+no_sample_of_fullforecasts*H_step_forecast_max; % Total forecasts
Index_full_forecast=1:H_step_forecast_max; % Index of full forecasts
Index_full_forecast_no_samples=repmat(Index_full_forecast,no_sample_of_fullforecasts);
Index_full_forecast_no_samples=Index_full_forecast_no_samples(1,:); % Samples with full forecasts and index

no_full_forecats=repmat(H_step_forecast_max,no_sample_of_fullforecasts); % Number of full forecasts
no_full_forecats=no_full_forecats(1,:);
no_nonfull_forecats=1:H_step_forecast_max-1; % Number of non_full forecasts
no_forecasts=[no_nonfull_forecats no_full_forecats]; % Number of forecasts (total=full+non_full)
no_forecasts_cumsum=cumsum(no_forecasts,2); % Cumsum Number of forecasts

Index_nonfull_forecast=zeros(Total_forecast-size(Index_full_forecast_no_samples,2),1)'; % Index of non_full forecasts
for i_H=1:H_step_forecast_max-1
    if i_H==1
  Index_nonfull_forecast(:,no_forecasts_cumsum(:,i_H))=1; 
    elseif i_H>1
 Index_nonfull_forecast(:,no_forecasts_cumsum(:,i_H-1)+1:no_forecasts_cumsum(:,i_H))=1:i_H; 
    end
end
Index_forecast=[Index_nonfull_forecast Index_full_forecast_no_samples]; % Index of forecasts

% Create forecast error matrices
FE_full_VAR=NaN(n,Total_forecast); % For VAR
FE_full_LP=NaN(n,Total_forecast);  % For LP

% Run Cross-validation forecast
for i_t=1:T_half
 
% Create traning and test sets
y_train=data(1:end-i_t,:);
y_test=data(end-i_t+1:end,:);

% Correct test set
if i_t<=H_step_forecast_max
    H_step_forecast=i_t;
    y_test_check=y_test;
else 
    H_step_forecast=H_step_forecast_max;
    y_test_check=y_test(1:H_step_forecast,:);
end

% VAR forecasting
[~,y_train_var_lag1,y_train_var_ad]=lag_variable(y_train,p);
B_VAR_train=(y_train_var_lag1'*y_train_var_lag1)\(y_train_var_lag1'*y_train_var_ad);

% Iterated forecasting
X_T_var=[1 y_train_var_ad(end,:) y_train_var_lag1(end,2:end-n)]';
for jj = 1:H_step_forecast    
    yhat_var(jj,:) = B_VAR_train(1,:)' + B_VAR_train(2:end,:)'*X_T_var(2:end);
    
    % Reconstructing the lag matrix with the new forecasts
    X_T_var = [1; yhat_var(jj,:)'; X_T_var(2:end-n)];
   
    % Square forecast error
    FE_VAR(jj,:)=(yhat_var(jj,:)-y_test_check(jj,:)).^2;   
end

% Store square forecast error
if i_t==1
FE_full_VAR(:,no_forecasts_cumsum(:,i_t))=FE_VAR';
elseif i_t>1
FE_full_VAR(:,no_forecasts_cumsum(:,i_t-1)+1:no_forecasts_cumsum(:,i_t))=FE_VAR';
end
FE_VAR=[];

% LP forecasting
B_LP_train=(y_train_var_lag1'*y_train_var_lag1)\(y_train_var_lag1'*y_train_var_ad);

% Sequence of coefficient matrices
for h=2:H_step_forecast
    Xh_train_projSet=y_train_var_lag1(1:end-h+1,:);
    Xh_train=y_train_var_ad(h:end,:);
    Bh_LP_train(:,:,h-1)=(Xh_train_projSet'*Xh_train_projSet)\(Xh_train_projSet'*Xh_train);
end
X_T_lp=[1 y_train_var_ad(end,:) y_train_var_lag1(end,2:end-n)]';

% Direct forecasting
for jj = 1:H_step_forecast   
     if jj==1
     yhat_lp(jj,:) = B_LP_train(1,:)' + B_LP_train(2:end,:)'*X_T_lp(2:end);
     elseif jj>1
     yhat_lp(jj,:) = Bh_LP_train(1,:,jj-1)' + Bh_LP_train(2:end,:,jj-1)'*X_T_lp(2:end);
     end

     %Square forecast error
     FE_LP(jj,:)=(yhat_lp(jj,:)-y_test_check(jj,:)).^2;   
end

% Store square forecast error
if i_t==1
FE_full_LP(:,no_forecasts_cumsum(:,i_t))=FE_LP';
elseif i_t>1
FE_full_LP(:,no_forecasts_cumsum(:,i_t-1)+1:no_forecasts_cumsum(:,i_t))=FE_LP';
end
FE_LP=[];

end % End cross-validation forecasting

% Compute Root Mean Square Error
% Create matrix for MSE
MSE_each_hor_index=1:H_step_forecast_max;
MSE_each_hor=zeros(n,H_step_forecast_max);
MSE_each_hor_lp=zeros(n,H_step_forecast_max);

%Compute MSE method 1: (MSE for each h-step ahead forecast horizon only)
for i_hor=1:H_step_forecast_max
forcast_position=find(Index_forecast==i_hor);
FE_each_hor_VAR=FE_full_VAR(:,forcast_position);
FE_each_hor_LP=FE_full_LP(:,forcast_position);

MSE_each_hor_VAR(:,i_hor)=nanmean(FE_each_hor_VAR,2);
MSE_each_hor_LP(:,i_hor)=nanmean(FE_each_hor_LP,2);
end

%Compute MSE method 2 (MSE upto h-step ahead forecast horizon)
%for i_hor=1:H_step_forecast_max
%i_hor=2;
%i_hor_series=1:i_hor; 
%forcast_position = ismember(Index_forecast, i_hor_series);
%indexes = find(forcast_position);
%FE_each_hor=FE_full(:,indexes);
%FE_each_hor_lp=FE_full_lp(:,indexes);
%MSE_each_hor(:,i_hor)=nanmean(FE_each_hor,2);
%MSE_each_hor_lp(:,i_hor)=nanmean(FE_each_hor_lp,2);
%end

% Compute weight
k=1; % Exponent number (a scalar e.g 1 or 2,..)
inv_mse_var=(1./(MSE_each_hor_VAR)).^k;
inv_mse_lp=(1./(MSE_each_hor_LP)).^k;
sum_inv_mse=(1./(MSE_each_hor_VAR)).^k+(1./(MSE_each_hor_LP)).^k;

w_var=(inv_mse_var./sum_inv_mse);
w_lp=(inv_mse_lp./sum_inv_mse);

w_var_full=[0.5*ones(n,1) w_var]; % At h=0 and h=1, LP and VAR are the same so same weight
w_lp_full=[0.5*ones(n,1) w_lp];


end